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INTRODUCTION 


The  field  of  three-dimensional  {3D)  underwater  acoustics  has  been 
developed  in  an  attempt  to  describe  the  properties  and  characteristics  of 
sound  propagation  in  the  heterogeneous  ocean  environment  (Lynch  and 
Chiu,  1989).  At  present,  there  are  three  major  approaches  used  in  studying 
the  3D  underwater  sound  field.  One  of  these,  ray  theory,  is  premised  upon  a 
high-frequency  approximation  method.  This  method  gives  a  geometric- 
optics  solution  which  involves  ray-tracing  in  a  spatially  varying  sound  speed 
field.  The  Hamiltonian  Acoustic  Ray  Tracing  Program  for  the  Ocean 
(HARPO)  is  a  versatile  numerical  code  developed  by  Jones  et  al.  (1986)  for  the 
computation  of  3D  rays.  The  primary  deficiency  of  this  theory  is  that  it  cannot 
adequately  address  the  behavior  of  low-frequency  sound  due  to  the  neglect  of 
sound  dispersion  and  diffraction. 

The  second  approach  uses  a  parabolic  approximation  to  the  acoustic  wave 
equation,  which  was  introduced  by  Tappert  (1977).  A  3D  numerical  parabolic 
equation  (PE)  model  was  developed  by  Lee  et  al.  (1988)  using  an  implicit  finite 
difference  scheme.  Another  3D  PE  model  was  arrived  at  by  Baer  (1981) 
utilizing  a  split-step  Fourier  algorithm.  In  a  recent  study,  analytic  solutions  to 
the  3D  parabolic  equation  were  obtained  by  Seigmann  et  al.,  (1990).  These 
solutions  are  valuable  for  testing  the  accuracy  of  3D  numerical  models. 

A  third  approach  is  the  normal  mode  method.  Pierce  presented  a  three- 
dimensional  version  of  this  method  in  1965.  The  normal  mode  method  is 
based  upon  a  separable  solution  to  the  wave  equation.  Pierce  assumed  an 
adiabatic  acoustic  environment  which  leads  to  the  neglect  of  coupling 
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presented  by  Chiu  and  Ehret  (1990).  This  later  3D  normal  mode  model 
accounts  for  both  horizontal  refraction  and  model  coupling. 

Any  improvement  in  acoustic  modeling  must  be  quantified.  We  shall 
endeavor  to  test  the  accuracy  of  the  3-D  coupled  normal  mode  model  of  Chiu 
and  Ehret  (1990)  in  regard  to  horizontal  refraction. 

To  examine  the  accuracy  of  the  Chiu-Ehret  model,  we  compare  the  results 
of  their  model  with  those  from  analytic  solutions  to  the  parabolic  equation 
arrived  at  by  Seigmann,  et  al.  (1990).  The  results  of  the  Chiu-Ehret  model  and 
the  analytic  soluiions  of  the  parabolic  equation  are  compared  for  two  cases  of 
horizontal  variation  in  sound  speed.  As  a  means  of  comparing  the  results,  we 
examine  the  slow  \  ariations  of  the  complex  pressure  enx  elope  function  and 
transmission  loss. 

Our  simulated  acoustic  field  has  a  range  of  100  km  with  azimuth  from  30° 
to  180°.  Two  variations  in  the  horizontal  sound  speed  field  were  used.  In 
Case  I  the  sound  speed  varies  only  with  the  azimuth  angle  while  Case  II  varies 
both  with  the  azimuth  and  radial. 

The  normal  mode  theory  used  by  Chiu  and  Ehret  is  described  in  detail  in 
the  Appendix.  The  appendix  also  includes  a  description  of  the  PE 
approximation  and  the  anal)  tic  solutions  developed  by  Seigmann  et  al. 
(1990). 
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II.  ANALYSIS  AND  METHOD  OF  MODEL  COMPARISON 


The  three-dimensional  coupled  normal  mode  model  was  applied  to  two 
sound  speed  fields  for  which  exact  analytic  solutions  to  the  PE  approximation 
are  available.  These  analytic  PE  solutions  for  vertically  invariant  sound  speed 
were  developed  by  Seigmann  et  al.  (1990).  Their  development  was  based  on 
expressing  the  “time-independent"  acoustic  pressure  component  as 


0p£{r,d,z)  =  Up(r,9)sin(y^2)exp|- 


where 

k 

»  C„ 


/j  ”  ^"^2  H  0  ”  2/ •••) 


ikpr 

2 


2  {/.  n 

exp  i\^kQr-- 


Cq  =  reference  sound  speed 
/  =  acoustic  frequency 
H  =  ocean  depth. 

The  reduced  complex  envelope  function 


Up{r,e)  =  Ap{r,e)exp{iG{r,e)} 


lip  is  governed  by 


where 
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«{r)  =  ^  {f  =  {r,d,z)) 

C(r)  =  sound  speed  as  a  function  of  position  r. 


Substitution  and  separation  of  real  and  inaaginary  parts  gives  two  coupled 


equations  governing  and  0: 


dAp  I  dAp  30  Ap  0^0 
dr  JcQr^  do  do  2kQr^  dO^ 


1 

+ 

1 _ 

1 - 

_ 

\  ) 

21 

UoJ 

'j 

1 

2 


2  ae  _L__L^^ 

*:o  Sr  *  fide)  (fcorf  ' 


(2.2) 


(2.3) 


If  one  chooses  0(r,O)  =  a(r)0,  equation  (2.2)  may  be  solved  for  the  amplitude, 
i.e., 


Ap(r,e)  =  F 


1  ra(s) 


to'  s' 


ds 


(2.4) 


Choosing  a(r)  and  the  functional  dependence  F  results  in  a  reduced 
envelope  function  (and  ultimately  acoustic  pressure)  while  equation  (2.3)  can 
be  solved  for  the  index  of  refraction  n(r)  and  thus  the  sound  speed  field  C(f). 
Equation  (2.3)  is  for  a  wide  angle  PE  approximation.  This  wide  angle  form 
will  be  used  throughout.  To  judge  the  accuracy  of  the  Chiu-Ehret  model  in 
calculating  horizontal  refraction  we  compare  the  numerical  coupled  mode 
results  with  the  analytic  PE  results  for  two  sound  speed  fields  or  two 
functional  definitions  of  a(r). 
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A.  INDEX  OF  REFRACTION 
1.  Case  I 

For  Case  1  we  choose  a(r)  =  Oq/Tq  r,  or 

0  =  ocqIcq  rd  (2.5) 


where  Uq  is  an  arbitrary  constant.  From  equation  (2.4)  the  amplitude  is 

A^{r,e)  =  Fie-a^W). 

Let  the  function  F  be  multiplied  by  The  amplitude  has  the  final  form 

A^(r,e)  =  (iQb^ie-a^lnr)  (2.6) 


where  ^  is  an  arbitrary  constant  with  units  of  length.  The  square  of  the  index 
of  refraction  from  (2.5)  and  (2.6)  is 


=1 +2 


c  l 

2' 

- 

<•  N 

2‘ 

2  r 

l.i 

71  ' 

4 

1  +  - 

lx 

•  2 

2 

.  ^0. 

2 

.^0. 

4ao0  +  IttQ  + 


II 

L  ^-0  J 


.  (2.7) 


The  resultant  sound  speed  is  shown  in  Figure  1  for  =  0.005  and 
Cq  =  1500  m/s. 

2.  Case  II 

For  Case  II  we  choose 

0  =  ayor^e  (2.8) 


and  hence 

A^,{r,e)  =  pQk^^(e-aQkQr).  (2.9) 

The  square  of  the  index  of  refraction  is  now 
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The  resulting  sound  speed  is  shown  in  Figure  2  for  Uq  =4x10  and  Cg  =  1500 
m/s. 

B.  ENVELOPE  FUNCTION 

We  now  have  two  analytic  sound  speed  fields  and  the  associated  acoustic 
pressure  fields  derived  from  the  parabolic  equation  approximation. 
Considering  only  j  =  1,  Seigmann  et  al.,  (1990)  define  pressure  as 

1 


0pE{r,e,z)  = 


ynkQTj 


^  2 


ikQV 


ll 


1-1 

4 


(2.11) 


For  Case  I,  using  the  amplitude  and  phase  expressed  in  (2.6)  and  (2.5), 
respectively,  we  obtain 


0p£{r,O,z)  = 


VTrfcory 

xexp 


2  ^2  ‘  *^0'' 


-iVfli 
2  Uo 


^2 


^2 

1 

4^0/ 


i-lii 


(2.12) 


For  Case  II,  using  (2.8)  and  (2.9),  we  get 
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0PE{r,9,z)  = 


^sin(yi2) 


I 

fiif 

1-i 

4 

fn] 

2T 

l^oJ 

(2.13) 


In  terms  of  normal  modes,  pressure  can  be  expressed  as 


‘»NM(r,e.z)  =  4.XU„{r,eyK*"<’''®>*Z„(z:r,e)  (2.14) 

n 

th 

where  is  the  n  normal  mode  (eigenmode),  is  the  corresponding 
horizontal  wavenumber  (eigenvalue)  and  L7„  is  the  corresponding  slowly 
varying  modulation  envelope.  See  the  appendix  for  details. 

Since  we  are  considering  depth-inva»-’ant  sound  speed  fields,  an  analytic 
solution  is  available  for  the  eigenvalues  and  eigenmodes.  The  first  mode  is 


and  thus,  keeping  only  the  contribution  of  the  first  mode,  (2.14)  becomes 

^NM{r>0,z)  =  ^ ^ til {r,^)e  ^0*^1  ('■'^)‘^''sin(7iz).  (2.15) 

The  normal  mode  model  entails  a  numerical  calculation  of  the  slowly 
varying  envelope  function  Therefore,  the  quantification  of  model 

accuracy  can  be  achieved  by  comparing  the  numerical  to  those  derived 
from  the  analytic  PE  solutions.  The  analytic  envelope  function  U^{r,6)  for 
Case  I  is  derived  by  equating  (2.12)  to  (2.15).  The  resulting  expression  is 


8 


H 


Uiir,0)  =  J—j3oko(e-aolnr)expi 


kor 


1-1 

2 


2r 


1-1 
4 


II 


+  ao0 


(2.16) 


where 


kf  =k^-  Yi 
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M) 

With  the  index  of  refraction  n  given  by  (2.7),  the  range  integral  of  k-^  can  be 
expressed  as 


jk-^dr=kQrl 
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It  follows  that  the  analytic  envelope  function  U-^(r,6)  can  be  recast  as 


U,(r,e)  = 


ft)(:o(e-Oo  >">')<! 


(2.17) 


where 
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+  cxq  6  — ' 

3-, 

4 

1  +  1 

11 

M) 

-2 

4ftQ0  +  2o!Q  + 

f 

Kh) 

4 

1 

2 

. 
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For  Case  II,  the  wave  number  integral,  using  (2.10),  becomes 
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III.  MODEL  RUNS 

Our  computational  domain  is  in  a  cylindrical  coordinate  system.  The 
horizontal  range  is  from  1  km  to  100  km.  The  azimuth  range  is  30°  to  150°. 
The  depth  is  4000  m.  The  domain  is  the  same  for  both  cases. 

We  use  a  frequency  of  50  Hz  and  source  depth  of  1000  m.  The  reference 
sound  speed  is  1500  m/s.  We  consider  mode  1  only.  The  vertical  boundary 
conditions  for  the  problem  are  pressure  release  surface  and  rigid  bottom.  Two 
range  and  azimuth  dependent  sound  speed  fields  as  presented  in  the  last 
chapter  are  used.  For  Case  I  sound  speed  ranges  from  1481  m/s  to  1496  m/s, 
varying  only  in  azimuth  with  a  gradient  of  0.128  to  0.130  m/s/degree  (see 
Table  1).  For  Case  II  sound  speed  varies  in  both  range  and  azimuth  (see  Table 
1  again).  Here  sound  speed  goes  from  1437  m/s  to  1500  m/s  with  a  radial 
gradient  of  0.13  to  0.657  m/s/km  and  an  azimuth  gradient  of  0.004  to  0.431 
m/s/degree.  Both  selections  of  sound  speed  fields  contain  realistic  gradients 
as  observed  in  the  ocean. 


TABLE  1.  THE  AZIMUTH  AND  RADIAL  SOUND  SPEED  VARIATION 


Case 

Range 

r{km] 

Sound  Speed  [m/s] 

Azimuth  Variation 

Radial  Variation 

f 

e=30° 

0=90° 

0=150° 

0=30° 

0=90° 

0=150° 

0=30° 

0=90° 

0=150° 

I 

I 

1496 

1488 

1481 

0.1303 

0.1289 

0.1275 

0 

0 

0 

50 

1496 

1488 

1481 

0.1303 

0.1289 

0.1275 

0 

0 

0 

100 

1496 

1488 

1481 

0.1303 

0.1289 

0.1275 

0 

0 

0 

II 

1 

1500 

1500 

1499 

0.0044 

0.0043 

0.0043 

0.1316 

0.3944 

0.6482 

50 

1493 

1480 

1468 

0.2217 

0.2178 

0.2142 

0.1310 

0.3848 

0.6221 

100 

1487 

1462 

1437 

0.4265 

0.4127 

0.4001 

0.1301 

0.3763 

0.5992 
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Eigenmodes  and  eigenvalues  are  calculated  using  finite  difference 
approximations  to  (A.8)  in  the  appendix  with  a  40  m  vertical  grid  spacing  and 
a  matrix  eigenvalue  solver. 

To  obtain  the  numerical  envelope  function  a  first  order  differential 
equation  ((A.24)  in  the  appendix)  is  integrated.  The  radial  integration  step 
size  is  0.5  km.  The  61  integration  paths  were  separated  by  two  degrees. 
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IV.  RESULTS 


A.  CASE  I 

In  comparing  the  numerical  normal  mode  (NM)  and  analytic  parabolic 
equation  (PE)  results  we  focus  on  the  slowly  varying  modulating  pressure 
envelope.  The  removal  of  the  rapid  oscillations,  which  have  wave  lengths  of 
order  27r/kQ,  makes  the  displays  of  results  easier  to  interpret. 

Figures  3  and  4  show  the  amplitude  and  phase,  respectively,  of  the 
envelope  function  calculated  by  the  normal  mode  model  using  the  Case  I 
sound  speed  field.  Figures  5  and  6  are  the  amplitude  and  phase,  respectively, 
of  the  envelope  based  on  the  analytic  solution  to  the  parabolic  equation. 
Figure  7  shows  the  relative  difference,  defined  as 

Up£(  analytic  solution)  -  (numerical  solution) 

U|,£  (analytic  solution) 

The  percent  difference  is  everywhere  less  than  5%.  The  difference  in  NM  and 
PE  phase  is  shown  in  Figure  8. 

Transmission  loss  at  a  depth  of  1000  m  calculated  by  the  NM  model, 
TLjvjn^^,  is  shown  in  Figure  9.  The  magnitude  of  the  analytic  PE  pressure,  from 
(2.12),  is 

''’'■S'''’'”*’'''’ 


Substituting  (2.6)  in  (3.1),  we  get 


|P|  = 


2 

\nkff 


-aolnr)sin 
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The  analytic  PE  transmission  loss  TLp^  is  thus 

TLpE  =  -ZOlogiolf’l 

=  101ogio^^-201og|o|i3o^'o(^  -  -  201og 

Figure  10  displays  TLp^  at  a  depth  of  1000m.  Both  transmission  loss  results 
show  that  the  azimuth  variation  of  sound  speed  produces  low'er  loss  at  large 
azimuth  angle. 

The  difference  between  the  analytic  PE  transmission  loss  TLp^  and  the 
numerical  NM  transmission  loss  is  shown  in  Figure  11.  Since  the 

modal  transmission  loss  is  a  function  of  amplitude  alone  this  difference  has 
the  same  shape  as  the  relative  error  in  the  amplitude  of  the  envelope  function. 
The  difference  in  transmission  loss  is  everywhere  less  than  1  dB. 

B.  CASE  II 

In  the  second  case,  we  examined  a  sound  speed  field  that  varies  in  range 
and  azimuth,  closely  representing  an  eddy  or  ring  structure  in  the  real  ocean. 

Figures  12  and  13  show  the  amplitude  and  phase,  respectively,  of  the 
envelope  function  from  the  normal  mode  model.  Note  that  azimuth 
variation  of  the  amplitude  is  larger  than  range  variation.  Figures  14  and  15 
are  the  amplitude  and  phase  of  the  envelope  function  Upg  based  on  the 
analytic  solution  to  the  wide-angle  PE.  Figure  16  show's  the  relative 
difference  of  the  amplitude.  The  percent  difference  is  everywhere  less  than 
2%.  The  difference  betw'een  the  NM  and  PE  phases  is  show'n  in  Figure  17. 
The  difference  of  phase  is  everywhere  less  than  one  degree. 


.  TT 

sm - 2 

2H 
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Figure  11.  Difference  in  TL  in  dB  for  Case  I  with  a  -  0.005  and  =  0.1 
The  NM  and  PE  transmission  losses,  TLj^^  and  TLp^,  at  a  depth  of  1000  m 
are  shown  in  Figures  18  and  19,  respectively.  They  are  nearly  the  same  as  can 
be  seen  in  Figure  20  where  the  error  is  everywhere  less  than  0.2  dB. 


V.  DISCUSSION 


In  this  chapter  we  address  two  issues.  The  first  is  a  quantification  of 
horizontal  refraction.  The  second  is  an  assessment  of  model  accuracy. 

A.  HORIZONTAL  REFRACTION 

Critical  to  this  accuracy  test  is  the  use  of  sound  speed  fields  that  cause 
significant  horizontal  sound  refraction.  To  confirm  that  the  selected  sound 
speed  fields  do  cause  significant  horizontal  refraction,  we  compare  the 
corresponding  Nx2D  NM  solutions  (the  calculation  is  divided  into  N  vertical 
slices  and  a  2D  NM  model  is  used  in  each  radial  direction)  with  3D  NM 
solutions. 

The  envelope  functions  of  the  N  x  2D  solutions  have  been  calculated 
from  (A.25)  in  the  appendix.  The  amplitude  of  the  NM  envelope  function  is 
the  same  for  both  the  Nx2D  and  3D  solutions  in  Case  I  and  Case  II.  However, 
the  phases  are  different.  Figure  21  shows  the  phase  of  the  Nx2D  envelope 
function  for  Case  I.  The  phase  is  essentially  constant  at  -45°.  In  Case  II,  the 
phase  of  the  envelope  function  of  the  Nx2D  solution  is  also  constant  and  its 
value  is  also  -45°. 

Now  let  us  compare  the  phases  calculated  by  the  Nx2D  method  with  the 
phases  of  the  3D  solutions  for  Case  I  and  Case  II,  as  displayed  in  Figures  4  and 
13,  respectively.  Both  cases  show  a  13°-15°  phase  difference  over  a  range  of 
100  km,  implying  that  the  azimuthal  sound  speed  variation  used  is  large 
enough  to  induce  sound  propagation  out  of  the  vertical  plane. 
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Figure  21.  Phase  of  for  Nx2D  for  Case  I  with  a  =  0.005  and  P  =  0.1 

B.  ACCURACY 
1.  Phase 

In  this  section  we  discuss  the  differences  between  the  analytic  PE  and 
numerical  NM  solutions.  In  Case  I,  the  phase  difference  between  the  analytic 
PE  and  numerical  NM  solutions  is  less  than  2.8  degree  everywhere.  The  Case 
II  phase  differences  between  the  numerical  NM  and  analytic  PE  solutions,  as 
shown  in  Figure  17,  is  everywhere  less  than  1  degree.  A  quantification  of  the 
phase  differences  for  the  two  runs  is  given  in  Table  2. 


TABLE  2.  PHASE  DIFFERENCE 


Case 

PE 

Normal  Mode  NM 

Difference  1  PE-NM  1 

min 

max 

min 

max 

max 

Case  I 

-45.15 

-60.35 

-45.15 

-58.92 

2.8 

Case  II 

-45.00 

-59.04 

-45.00 

-58.94 

0.43 
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2.  Amplitude 

Results  are  shown  in  Table  3.  Envelope  amplitude  of  NM  and  PE 
generally  agree  well  with  the  largest  relative  difference  near  5%.  The  largest 
difference  occurs  in  a  region  that  has  the  largest  azimuth  gradient  in  sound 
speed. 


TABLE  3.  AMPLITUDE  DIFFERENCE 


Case 

Amplitude  of  NM 

Amplitude  of  PE 

Error  Value 
Difference 

Relative 
Difference  [%] 

min 

max 

min 

max 

max 

max 

Case  I 

0.798 

4.22 

0.817 

4.27 

0.038 

4.8 

Case  II 

0.851 

4.28 

0.841 

4.28 

0.076 

1.8 

3.  Transmission  Loss  (TL) 


The  difference  in  transmission  loss  between  the  analytic  PE  and 
numerical  NM  solutions  is  shown  in  Table  4.  The  difference  for  both  cases  is 
everywhere  less  than  half  a  decibel. 


TABLE  4.  TRANSMISSION  LOSS  DIFFERENCE 


Case 

TL  of  NM 

TL  of  PE 

Difference 

min 

max 

min 

max 

max 

Case  I 

28.84 

63.31 

28.85 

63.73 

0.43 

Case  II 

28.73 

62.75 

28.73 

62.85 

0.15 
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VI.  CONCLUSIONS 


In  this  work  we  have  tested  a  three-dimensional  numerical  acoustic 
coupled  mode  model.  We  have  sought  to  describe  its  accuracy  in  regard  to 
the  modeling  of  horizontal  refraction.  In  order  to  test  its  accuracy,  we 
compared  the  coupled  mode  model  results  with  two  different  analytic 
solutions  to  the  parabolic  wave  equation. 

We  have  concentrated  on  the  accuracy  of  the  envelope  function  and 
transmission  loss  calculation.  For  the  sound  speed  fields  described  in  Table 
1,  we  found  that  the  phase  error  of  the  slowly  varying  envelope  function  is 
lower  than  2  degrees  and  that  the  envelope  amplitude  agrees  to  within  5%  or 
better.  Our  test  cases  have  indicated  that  the  normal  mode  model  agrees 
closely  with  the  analytic  solutions. 

We  have  only  used  depth-invariant  sound  speed  fields  for  this  test.  As  a 
result,  the  accuracy  in  modeling  mode-mode  interactions  is  not  tested  here. 
Future  tests  should  be  directed  at  examining  the  accuracy  of  calculating 
mode-coupling  effects  using  depth  varying  sound  speed  fields. 
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APPENDIX 


where  C(r)  represents  the  speed  of  sound  propagation.  The  cylindrical 
coordinates  are  r  =  (r,  6,  z)  where  r  is  range,  9  is  the  azimuthal  angle 
measured  positive  counterclockwise,  -z  is  depth,  and  V  is  the  Laplacian 
operator  in  cylindrical  coordinates. 


1  a  1  a^  a^ 


(A.2) 


Inclusion  of  an  acoustic  source  modifies  Equation  (A. I)  to 

=  (A3) 

C  ot" 

For  an  acoustic  point  source  located  at  Fq  with  time-harmonic  circular 
frequency  w,  Equation  {A.3)  becomes 

=  -An5{f  - (A.4) 

C  dr 

If  we  let 
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p{r,t)=0(f)e-‘'"‘  (A.5) 

where  <P(r)  is  tiiiie-independent,  then  a  Helmholtz  equation  results: 

V^0{f)  +  k^{f)0{f)  = -47i5{r -fo)  (A-6) 

VO 

where  the  acoustic  wavenumber  is  fc(f)  =  .  For  r  the  Helmholtz 

c(r)  ° 

Equation  is 

V^0(f)  +  k^(f)0(f)  =  0.  (A.6a) 

b.  Coupled  Normal-  Mode  Solution 

Adiabatic  normal  mode  theory  was  applied  to  a  range-dependent 
medium  by  Pierce  [1965]  and  Milder  [1969].  The  pressure  field  was  expressed 
as  a  linear  combination  of  the  local  modes  whose  coefficients  were  obtained 
from  a  system  of  decoupled  ordinary  differential  equations.  The  physics  of 
coupled  normal  mode  theory  includes  non-adiabaticity  and  results  in  a 
coupled  system  of  differential  equations.  The  theory  used  in  the 
development  of  the  Chiu-Ehret  3D  coupled  normal  mode  model  follows. 

The  first  step  in  the  development  is  the  expansion  of  the  time- 
independent  pressure  component  in  terms  of  the  local  normal  modes,  i.e., 

<!>{r)  =  ^'LPn(’-,e)Z„{z:r,e)  (A.7) 

n 

where 

is  the  n^^  local  normal  modes  at  point  (r,9). 
is  the  n**’  mode  amplitude  function  at  point  (r,0). 

c.  The  Local  Normal  Modes  Z„ 

The  local  normal  modes  obey  the  following  equation: 
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9z 


Z„=0 


(A.8) 


th 


where  is  the  horizontal  wavenumber  associated  with  the  n  mode  at  each 
horizontal  location. 

The  idealized  upper  and  lower  boundary  conditions  are 

Z„(2  =  O;r,0)  =  O.  (A.9) 


— Z„(2  =  -H;r,0)  =  O. 
02 


(A.10) 


where  H  is  the  ocean  depth.  The  local  normal  modes  can  be  normalized 
according  to  the  orthonormal  condition 


H 


~  ^mn- 


(A.ll) 


d.  The  Mode  Amplitude  Function 

To  obtain  the  governing  equation  for  substitution  of  (A.7)  in 
(A. 6a)  is  necessary.  The  use  of  (A.8)  and  the  farfield  approximation 


1,2  11 

-^4^2 


mm  tn  m  m 


(A.12) 


following  the  substitution  results  in  the  following  equation; 


m 


,  1 

J  dPm 

dr^ 

r2  302  j 

■"lar  9r 

f d^z„  .  1  1 

+  k^  P 
^  m- 

t  9r2  r2 

902  J 

•  _  _  _ 


(A.13) 


The  next  step  is  to  multiply  (A.13)  by  Z^  and  then  integrate  over  z,  i.e.. 


32 


m 


|z»z, 


^  ^  ^  ■  1.2  p 

"  3r^  de^ 


dz  +  2  J  Z„ 


ap^  dZ„  ^  1  dPrr,  az^ 


+ 


dr  dr  dO  dO 


dz 


( d^Zm  1  a^z 


m 


dr^  r^  d9^ 


\dz 


=  0. 


(A.14) 


A  subsequent  application  of  the  orthonormality  condition  (A. 11)  on  (A.14) 
gives 


'  d^  ,2  1  a^  ^ 

dr^  r^  dO^ 


=-I 


m 


f  r 

2  jz„ 


az 


^dz 

dr  j 


2 


—  +  ±l  fz  ^ 
ar  ae 


ae 


+ 


Pmj^n 


( d^z^  1  a^z 


+ 


m 


ar^  d9^ 


dz 


Defining  coupling  coefficients  as 


J  J 
(A.15) 


y„„  =  2\z„^dz 

B  =-fz  -^-^-dz 

Pmn  ^  J 


Pmn  ~  j 


d9 

f  d^Zrr,  1  a^Z 


ar' 


d9^ 


dz, 


the  coupled  system  (A.15)  can  be  recast  as 


2  -^,32 

dr  r  d9 


p»  =  -Z 


m  i- 


i-  R  lA 

^"•"dr  '^""’rde 


+  B. 


mn 


(A.16) 


e.  The  Envelope  Functions 

Following  the  work  of  Chiu  and  Ehret  (1990),  the  mode 
amplitude  function  can  be  separated  into  a  slowly  varying  envelope 
function  and  a  rapidly  varying  component  i.e.. 
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0n=jlK{r,O)dr. 

A  substitution  of  (A.  17)  into  (A.  16)  yields 


(A17) 

(A18) 


^  ,2  1 
+  kt  + 


dr 


2  ”  r2  302 


=-S  +  uj* 


which  can  be  expanded  as 


'"ii 

dr^  dd^ 


+i2k  -^+J/^+X^^-Ar^YLl[^A 


^  dr  1  dr  dd^  r^\dd)\  r^  dd  apj  ” 


^2' 


(A19) 


=  -X 


m 


ymn  ^  +  Pmn  f  ^yrnnk-m  +  iPmn  — ^ 


mn 


dr  ’^'"'rdd)  30 

Dividing  (A.20)  by  gives  the  governing  equation  for  the  envelope: 


(A.20) 


Dr 


xdr"^  r^  dd^  j 


Ufi  + 


|;  +  £„ju„  +  F„i^U„  =  £G„„VU„H-XH„„U„ 


m 


m 


(A.  21) 


where 


Dr 


-1 

2L, 


2kr 


dkfi  ^  1  d^0n  ^  i  fd0n^ 


dr  302  ^2  00 


j_r^' 

k„r{  39  , 


^mn 


yymnf  +  P,.ney^*"~*'''> 
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V 


^Snn  im  =  n) 

ZKn 

and  r  and  6  are  the  unit  directional  vectors  in  r  and  6,  respectively. 

/.  Transmission  Loss  TL 

The  mode  amplitude  function  is  related  to  the  envelope 
function  and  the  phase  by  (A.17).  The  local  normal  modes  are 
solutions  to  the  eigenvalue-eigenfunction  problem  (A.8)-(A.ll).  Once  and 
Z„  are  computed,  pressure  can  be  calculated  from  the  product  of  and  Z^. 

At  each  point  (r,  6,  z),  0  is 

<P“(r,0,2)  =  [Re(«P)f  +[lm(<P)f  (A.22) 

where 

Re[<t(r,e,z)]  =  £[Re(P„)].Z„(2) 
m 

and 

m 


In  (A. 22),  Re  and  Im  are  used  to  denote  the  real  and  imaginary 
parts,  respectively.  Transmission  loss  can  be  computed  as 

TL  =  -10  log,o<P^fr,e,z).  (A.23) 
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2.  Numerical  Solution 

a.  Numerical  Model 

To  obtain  a  normal  mode  solution  we  must  solve  the  equations 
for  the  eigenvalues  and  eigenmodes  (A.8)-(A.ll)  and  the  equation  for  the 
modulation  envelope  (A. 21).  These  equations  can  be  solved  numerically. 

b.  The  Eigenvalues  and  the  Eigenmodes 

Equation  (A. 8)  is  approximated  using  central  finite  differencing. 
This  approximation  casts  (A.8)-(A.10)  into  a  matrix  algebra  problem.  The 
eigenvalues  and  the  eigenmodes  of  which  can  be  determined  using  the 
iterative  QR  method  (Acton,  1970). 

c.  The  Iterative  Method 

To  obtain  the  envelope  function  L/,,,  we  need  to  solve  equation 
(A. 21)  which  can  be  rearranged  to  form  a  set  of  first  order  partial  differential 
equations  (PDE)  with  smaller  terms  put  on  the  right-hand  side.  We  can  solve 
for  in  an  iterative  fashion.  The  iterative  equation  is 


—  +  E, 

lar  ' 


1  a 


U‘  -  y  -  D„ 

II  Li  III  ”  f.  II  ” 


III 


dr^ 


1  a^ 


U 


i-l 


where  u'„  is  the  solution  at  the  i‘'‘  iteration  for  mode  n.  This  set  of  first  order 
differential  equations  is  integrated  using  a  Runge-Kutta  method  of  order  5 
and  6  (Acton,  1970). 


3.  The  Nx  2D  Method 

In  the  N  ^  2D  method,  the  sound  channel  is  divided  into  N  vertical 
slices  and  the  two-dimensional  normal  mode  solution  is  solved  for  each  slice. 
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It  combines  the  results  in  each  vertical  plane  to  construct  a  3D  field.  This 
method  provides  an  approximation  to  a  3D  solution  for  a  3D  sound  field. 

Thus,  in  each  slice  the  envelope  function  U is  independent  of 
azimuthal  variation  effects.  With  the  azimuthal  terms  ignored,  the  governing 
equation  becomes 


dr 


1-1 


K  ■  I  ■  -D,, ^  I  G„,„VU';' 

m  m 


(A25) 


where 


D.,  = 


-1 

2E 
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'}}in 


1  dk, 


H 
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2k,j  dr 

ZK,i 

d  . 

—  r 

dr 

[jn  *  n) 

_L_  R 

2k„ 

{m  =  n) 

B.  THREE-DIMENSIONAL  PARABOLIC  WAVE  EQUATION 
1.  Parabolic  Wave  Equation  (PE) 

The  following  work  was  developed  by  Tappert  (1977)  and  Lee  (1988). 
In  cylindrical  coordinates,  the  Helmholtz  Equation  (A. 6a)  becomes 


d^0  d0  ^  1  d^0 
dr^  r  dr  d^ 


+  kln^0{r,e,z)  =0 


(A.26a) 


where 
A'(,  =  n>/Co 


37 


Cq  is  a  reference  sound  speed 
xv  is  the  harmonic  source  frequency 

n  is  the  index  of  refraction 
n  =  n(r,b,  z)  =  Cq  /C(r,6,  z). 

The  PE  approximation  assumes 

0{r.e,z)  =  u(r,e,z)tir) 


where  v{r)  is  rapidly  varying  and  u(r,  0,  2)  is  a  modulation. 
By  substituting  (A. 26b)  into  (A. 26a),  one  obtains 


'd\i 

fi  2at0 

,  r  V  dr  ^ 

dll 

dr 


+ 


1  d^u 
d6^ 


d^u 


+  k^n^u 


■ 

d^v 

1  du 

V  + 

4-  _  -1— 

rdr^ 

u 


If  V  is  determined  by 


d^v 

dr^ 


+ 


'i-dZ’2  r, 
-—■^k^v  =  0 
r  dr 


then  u  must  satisfy 


d^ii 

dr^ 


1  ^  2 
r  V  dr 


dit  1 


d\i 


dr 


d9^ 


dz 


j  +  /cq  •  l)i(  =  0. 


The  solution  to  (A. 27)  is 


v(r)  =  Hj(/cor)  ~ 


1 


f 

2 

2 

“ 

r  >  - 

^nk^r^ 

exp 

i 

‘o’- -7 

^  *  J  - 

where  Hq  is  a  Hankel  function  of  zeroth  order. 

By  substituting  (A. 29)  into  (A. 28),  one  gets 


iA.26b) 


0.  (A.26C) 


(A.27) 


(A.28) 


(A.29) 


38 


3^  d  3"  1  3^  ,2/  T'v 


u  =  0 


(A.30) 


which  can  be  factored  into 


fd 


U  =  0 


(A.31) 


where 


Q  = 


J_ii 

fco  ^2^ 


+ 


_1 _ ^ 

)t^r2  33“ 


2 

2 


Q  is  called  the  square  root  operator.  If  one  considers  only  outgoing  waves,  the 
equation  to  be  solved  is 


—  +  ikQ  -  ikoQ 
dr  j 


M  =  0. 


(A.32) 


Let 

1 

Q  =  [l  +  X  +  y]I  (A.33) 

so  that 


y  = 


1  3^ 

kl?  33^ 


and  then  expand  Q  in  a  Taylor  series; 

Q  =  [i  +  x  +  y]i  =  i+-x--x^+-y+.... 


(A.34) 
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Inclusion  of  the  first  four  terms  in  the  expansion  leads  to  a  "wide  angle" 
version  of  the  PE  approximation.  Substituting  (A. 34)  into  (A. 32)  gives 


dr  2  '  '  dz'^ 


{A.35) 


2.  The  Envelope  Function  of  the  Parabolic  Equation 

Here  we  present  the  analytical  results  of  Seigmann  et  al.,  (1990).  Their 
development  uses  the  pressure  release  surface  and  rigid  bottom  boundary 
conditions,  i.e.. 


u(r,0,z  =  0)  =  0 

^(r,0,2=  H)  =0. 
dz 


(A.36) 


A  solution  to  equation  (A.35)  for  a  depth-independent  sound  speed  field  and 
subject  to  the  boundary  conditions  stated  above  can  take  the  following  form: 


i/(r,0,2)  =  L/j;(r,0)sin(y^2)exp 


iV  II 
2  .  ^'0 , 


(A.37) 


where 


yj 


n 

H' 


H  is  the  ocean  depth,  and  U ^  is  the  envelope  function  of  the  parabolic 
solution. 

By  substituting  (A.37)  into  (A.35),  one  obtains  the  equation  governing 
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SU,,(r.9)  _  ll'  ifZi.]\_l_iLL„o, 

ar  2  I"  4!"  ■')  -zitoj  (V)^ae2 


Letting 


U^(r,0)  = 


(A.38) 


(A. 39) 


where  both  A  (r, 6)  and  0(r,0)  are  real  quantities,  (A.38)  becomes 


do  'M/;  ( 2  ,  I  yj  1  f  2  I  yj 

dr  dr  2  '  '  2  Kq  4^  ^  h 


+  2/^— +  jA„— ^  -  —  ,4 


i 

‘a^A^ 

2A'or 

L  aa" 

ae  ae  I  ae 


{/\.40) 


A  separation  of  the  real  and  imaginary  parts  of  (A. 40)  yields 

1  ^y^p  a©  A;  a^©  ^ 

-  +  - ~ - +  - r - TT  =  0 

dr  k^r^  dd  d9  2k^r  d^ 


(A.41) 


(2  h2  1  7,  2  a©  1  fa©l^  1  1  ^  A, 

[  2[A-oJ  J  4^  2[A:oJ  k^  dr  (A-pr)n^0j  [k^rf  y^p  dO^ 

(A.42) 

Given  and  ©,  a  pressure  field  is  determined,  and  the  index  of  refraction  can 
be  calculated  from  (A.42).  The  resultant  sound  speed  field  can  then  be  used 
in  the  normal  mode  model. 
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